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ABSTRACT 

I review the multiphase coohng flow equations that reduce to a relatively simple form 
for a wide class of self-similar density distributions described by a single parameter, k. 
It is shown that steady-state cooling flows are not consistent with all possible emissiv- 
ity profiles which can therefore be used as a test of the theory. In combination, they 
provide strong constraints on the temperature profile and mass distribution within 
the cooling radius. The model is applied to ROSAT HRI data for 3 Abell clusters. 
At one extreme (fc ~ 1) these show evidence for cores in the mass distribution of size 
110-140 /i^Q^kpc and have temperatures which decline towards the fiow centre. At the 
other (fc i-^ oo), the mass density and gas temperature both rise sharply towards the 
flow centre. The former are more consistent with measured temperatures, which sug- 
gests that the density distribution in the intracluster medium contains the minimum 
possible mixture of low-density components. 
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1 INTRODUCTION 

The multiphase nature of the intracluster medium (ICM) in 
cooling flows was demonstrated a decade ago when deprojec- 
tions of X-ray surface-brightness proflles showed that mass 
cools and is deposited from the flow in a distributed man- 
ner, M 2; r (e.g. Thomas, Fabian & Nulsen 1987). However, 
the complexity of the theory and lack of data with high spa- 
tial resolution means that the single-phase approximation is 
still widely adopted. In this paper I show how the use of self- 
similar density distributions can lead to a relatively simple 
form for the multiphase cooling flow equations whilst span- 
ning the whole range of expected behaviours in the more 
general case. In combination with the emissivity profile (for 
a spherically-symmetric flow), the equations can be solved 
to yield the gas temperature and mass proflles within the 
cooling flow. 

In Section ^, I introduce the concept of the volume frac- 
tion to describe the distribution of density phases in the 
ICM. The self-similar forms of the volume fraction are deter- 
mined and shown to span all reasonable behaviours for the 
more general situation. The steady-state, self-similar form of 
the cooling flow equations are derived in Section |^ and it is 
shown how these can be solved in the spherically-symmetric 
case if the emissivity proflle is known. The behaviour of the 
solutions is examined in detail: not every emissivity profile 
is consistent with these models. In Section ^, the theory is 
applied to ROSAT HRI data for three Abell clusters with 
some success. Finally, the results are summarised and fur- 
ther discussed in Section ^. 



2 THE FORM OF THE DENSITY 
DISTRIBUTION 

The theory of multiphase cooling flows was set out by Nulsen 
(1986). He introduced the concept of the volume fraction, 
f{p,r,t), such that f dp is the fractional volume occupied 
by phases with densities in the range p to p + dp. It is as- 
sumed that the phases comove but are thermally isolated 
from one another — these can be regarded as empirical facts 
as otherwise the flow would rapidly evolve to a single-phase 
state. 

Given these assumptions then it is relatively straight- 
forward to derive an equation for the evolution of the volume 
fraction (see Appendix). Writing 



(2 — q) (4-a)/(2-a) , ,n 



PO 



where 
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and p — p{r,t), then 
£ + (3-a)^-^ 
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V.u = 



(1) 



(2) 
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where the dot refers to a covariant derivate, following the 
flow. 

In general <; is a complicated function of position and 
time. However, we can look for solutions in which g has 
a constant functional form, g — g{w). Only the flrst term 
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in equation |^ depends upon w. Hence we require that g = 
wdg/dw oc g. There are two kinds of solution; 

(i) (?oo oc exp(~ui). This is the most extended distribution 
which is convectively stable (it gives P oc pg). It includes 
phases of arbitrarily low density. 

(ii) oc (1 — w) < w < 1; k > 1. These solutions 
possess a minimum density, p > po- k — 1 is the least ex- 
tended, consisting solely of the power-law cooling tail. As 
fc oo the solutions resemble Qoo- 

For other forms of g one must resort to numerical inte- 
gration to follow their evolution. Thomas (1988b) looked at 
the steady-state evolution of a range of distributions with 
a sharp cut-off at low densities and reached the following 
conclusions: 

• All distributions develop a high-density tail, / ^ 

^ as they cool. 

• Sufficiently narrow distributions resemble the pure 
power-law gi by the time they begin to be deposited. 

gi and goo bound all reasonable solutions of the cooling 
flow equations, be they self-similar in form or not. In prin- 
ciple, k could be less than unity (as i— » the flow tends 
to the homogeneous case) but it is difficult to see how such 
distributions could arise. If the large density variations in- 
ferred within the cooling radius result from amplification of 
an initially much narrower distribution, then values of k of 
order unity are to be preferred. Without a plausible forma- 
tion history for the ICM, however, it is better to leave k as 
a free parameter to be determined empirically. I will argue 
below that low values of k fit the observations of cooling 
flows in Abell clusters better than high values. 



3 RECONSTRUCTION OF CLUSTER MASS 
PROFILES 

The self-similar density distributions derived in Section ^ 
lead to particularly simple forms of the steady-state cooling 
flow equations. I derive these below and then show in the 
next section how they can be combined with the emissivity 
profile to provide strong constraints of the mass distribution 
within the cooling radius. 

3.1 The equations 



A2 



A3 



Substituting the functional form of gk into equations 
and 1^ we see that the self-similar forms of the cooling flow 
equations are: 



po 
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To these may be added the equation of hydrostatic support 
VP 



where <1> is the gravitational potential. I assume here that 
the inflow is subsonic — this turns out to be a good approx- 
imation in all multiphase cooling flow models. 

In a steady-state and spherical symmetry the above 
equations reduce to 

dlnP 

-r, = -2S, 8 

dmr 

dlnp 



dlnr 7 



and 
dlnti 
dlnr 

where 

^ _ GM pmn 
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2r ksT 

is the ratio of the virial to the thermal temperatures, and 
7-1 rigA(ro) r _ 1 dlnM 



(2 



7 



u fc dlnr 



(12) 



is the ratio of the inflow time to the constant-pressure cool- 
ing time. 

Although there appear to be three equations here, the 
dimensionless ratios S and r are the only important vari- 
ables. The third equation merely acts as a scaling (for fixed 
r, p^~" oc p^~°'u). Hence the physics can be captured in just 
two equations: 



dluE 

d In r 

and 
dlnr 
d In r 



X - 



= 3 



7-1 
1 + 2- E 



-[(3- 
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-7(l-a)]E- 



(13) 



(14) 



where x = dluM/dlnr. The goo equations can be recovered 
by letting fc oo and using fcr in place of r as the second 
dimensionless variable. 

The usual way of proceeding when solving the single- 
phase equations is to pick a functional form for the mass 
profile, Xi S'lid mass deposition rate, 13, then to solve for 
E and r. From this one can generate an emissivity profile, 
^(r), for comparison with the data. Here I adopt a different 
approach: because (3 is fixed (for a particular value of k) in 
the multiphase models, one can specify ^(r) and determine 
the mass profile. 

Suppose that /3flt = — (l/6)dln^/dlnr is known. Then 



dlnr 
d In r 



= 3 



, (3 — a) — 7(1 — a) 
2 — a + a'y 

27 

(2-a)(2-a-ha7) 



(15) 



0, 



(7) 



Furthermore this is an eigenvalue problem: requiring that 
the solution extend to r = fixes the outer boundary con- 
dition. Hence we can solve for r, E and x- 



3.2 The behaviour of solutions 

We can get a good idea of the behaviour of the solutions 
to equation |l^ by looking at the case of constant /3fif Im- 
posing the physical constraints E > (i.e. a non-negative 
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temperature) and x ^ (i.e. mass constant or increasing 
with radius) restricts /3fit to lie in the range 

80 + 21fc 



2(5 + 3 fc) 



120 + 36fc 



(16) 



(in this expression and henceforth I set 7 = 5/3 and a = 0.5 
rather than including them explicitly). Thus steep emissiv- 
ity profiles, /3fit ^ 0.65, are incompatible will all steady-state 
cooling flow models (larger values can occur, however, out- 
side the cooling radius). In addition, the inner value of 
can be used to constrain k: flat cores are inconsistent with 
small values of k and a central decrease in emissivity would 
be incompatible with all models. If we assume that the virial 
temperature drops (i.e. E 1— > 0) within the cluster core, then 
the central value of Pat provides a measure of k. 

The numerical solutions of equation |l3 in the case of 
non-constant Pat are illustrated in Figure nrHere I have set 
Pfit ~ 0.6 + 0.1 logjQ r and k = 1. Panel (a) shows the speed 
at which solutions diverge away from the desired one (i.e. 
the one which extends to r = 0) as r decreases. This ensures 
that the solutions are insensitive to the inner boundary con- 
dition and are stable to small variations in the value of Pat 
in the inner bin (which is most likely to be affected by the 
point-spread function, uncertain absorption correction, etc. 
). The difference in the two values of r at r = 10~^ in the 
flgure is less than one percent. Panel (b) shows that the so- 
lutions are able to cope with quite large variations in Pat'- 
in this case 100 fluctuations drawn from a uniform distri- 
bution of extent ±0.1 have been added to each decade in 
r. Despite the fact that Pat now varies outside both ends 
of the physical range described in equation it remains 
true to the smooth solution. In this paper we will flt smooth 
functions to observed emissivity profiles before attempting 
to solve the cooling fiow equations. 



4 APPLICATION TO ABELL CLUSTERS 

I will now apply the above theory to three cooling fiow clus- 
ters: A 85, A 745 and A 2025. The ROSAT HRI data were 
kindly supplied to me by Clovis Peres in the form of depro- 
jected density and temperature profiles. This has the advan- 
tage that conversion of the counts into emissivity, including 
correction for the spectral response and absorption, has al- 
ready been included. Provided that the cooling fiow solu- 
tions resemble the deprojected one, then these corrections 
will hold good. 

I will consider each of the clusters in turn. 



4.1 A 85 

The emissivity profile is shown in Figure ^. There are 23 
annular bins, each 12 arcsec wide, of which the inner 9 have 
cooling times less than 1.33 x 10^''yr (throughout this paper 
I take Ho = 50kms"^Mpc"^). The profile is well-fit by a 
simple broken power-law fit. 



0.12 Mpc 



+ 



0.12 Mpc 



(17) 




0.05 
r/Mpc 

Figure 2. The emissivity profile for A85 in 12 arcsec bins. The 
soUd fine shows a broken power-law fit as described in the text. 



will require an inner core in the mass distribution. This is 
illustrated in Figure |^. Note that E drops very close to zero 
at r = 10 kpc (it tends to a small constant value within this 
radius). The gravitational density profile (i.e. that of the to- 
tal mass, not just the gas) is well-fit at radii greater than 
20 kpc by a King model. 



0.11 Mpc 



(18) 



The asymptotic slope of ^ as fc i~> is very close to the min- 
imum permitted for fc = 1 which suggests that the solution 



Within this radius the density is poorly constrained. Al- 
though it appears to rise, only a small change in the slope 
of ^ would cause it to level off or even fall — all we know for 
sure is that the virial temperature becomes very small. Note 
also that there is only one bin within 20 kpc and this one is 
most likely to be affected by smoothing by the point-spread 
function, uncertain correction for excess absorption, etc. 

The temperature of the gas is approximately constant 
outside the core radius, but drops by a factor of five in to 
10 kpc. A temperature decline in the centre of clusters is 
typical of cooling flows observed by ASCA. 

Note that the slope of the mass-deposition proflle, r, 
is close to unity within the cooling radius, Tcooi ~ 150 kpc. 
This radius is not a special one for our solutions as we have 
assumed the cooling flow solution holds everywhere. For this 
reason the asymptotic slope of the gravitational density pro- 
flle at larger radii should be taken with a pinch of salt. 

The corresponding solution for A: = cxa is shown in Fig- 
ure ^. There in no core in the gravitational mass profile in 
this case, with pdark rising as all the way into 10 kpc. 
This is refiected in the temperature profile, however, which 
also rises by a factor of 5 between 150 and 10 kpc, unlike 
any observed cluster. 



4.2 A 745 

The emissivity profile, shown in Figure ^ has 19 bins of 
width 8 arcsec, with 9-12 bins interior to the cooling radius 
(''cool ~ 180-240 kpc). The first thing to note is that ^ fiat- 
tens considerably in the innermost bin. This is incompatible 
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-4 -3 -2 -1 -4 -3 -2 -1 

log, or log,or 

Figure 1. Solutions to the multipliase cooling flow equation, equation for simple test profiles: (a) /Jgt = 0.6 + 0.1 log^Q r, (b) 
/3fit = 0.6 + 0.1 logj^Q r ± 0.1. The curves which decline and increase with radius represent solutions for r and S, respectively. 




log^Q(r/Mpc) r/Mpc 

Figure 3. The A; = 1 cooling flow solution for A85; (a) E (dashed line) and t (solid line), (b) the density of the gravitating matter. 



in 




log.|g(r/l\/lpc) r/Mpc 
Figure 4. The A: = oo cooling flow solution for A85; (a) S (dashed line) and r (solid line), (b) the density of the gravitating matter. 
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0.01 



0.05 0.1 
logioCr/Mpc) 



Figure 5. The emissivity profile for A745 in 8 arcsec bins. Tlie 
solid line shows an analytic fit to the data and the dashed line 
shows the emissivity profile corresponding to a King-law density 
profile, as described in the text. 




0.05 0.1 

logioCr/Mpc) 



Figure 6. The emissivity profile for A2029 in 12 arcsec bins. The 
solid line shows an analytic fit to the data and the dashed line 
shows the emissivity profile corresponding to an analytic density 
profile, as described in the text. 



with density distributions with small values of k. In fact the 
fit shown by the solid line, 



1 + 



31 kpc 



+ 



89 kpc 



(19) 



admits solutions which extend to r = only for fe = oo. 
Unfortunately, just as for A 85, large values of k give steeply 
rising dark-matter density and temperature profiles at small 
radii. 

If we ignore the inner bin, however, then it is possible to 
find solutions for all values of k. The dotted line in Figure ^ 
shows the emissivity profile which corresponds to a mass 
density 



Pgra 



1 + 



0.11 Mpc 



(20) 



and k = 1. 
4.3 A 2029 

For A 2029 there are 7-9 bins of width 12 arcsec within the 
cooling radius, rcooiaj'P'"oa;170-210 kpc. The emissivity pro- 
file 

\ 0.0'± ^ 

r \ ( r 



0.13 Mpc 



+ 



0.13 Mpc 



(21) 



shown in Figure ^ is again too shallow in the central bin, but 
the inconsistency is this time so slight that i t st rengthens 
the case for fe = 1 (as discussed in Section 3.2 the slope 



of the emissivity profile within the cluster core provides a 
measure of k). Figure ^ shows that E drops to a value 
close to zero within 30 kpc. This indicates that the virial 
temperature has sunk well below the gas temperature (i.e. 
the flow is isobaric within this radius). This is reflected in 
Figure ^ which shows the corresponding density profile. 
The latter is well-fit at radii greater than 80 kpc by a mass- 
density 



0.14 Mpc 



(22) 



but declines rapidly within this radius. Only a small change 
in the emissivity profile is required to generate solutions in 
which the gravitational mass has a constant core density, as 
indicated by the dotted line in Figure ^ which reproduces 
the above density profile exactly. 



5 CONCLUSIONS AND DISCUSSION 

• I have rederived the cooling flow equations for self- 
similar density distributions. Two of these in particular are 
expected to bound the behaviour of all possible flows. 

• The steady-state cooling flow equations are not compat- 
ible will all conceivable emissivity profiles. Thus they can be 
used as a test of the theory. 

• The solutions provide bounds on Af (r) within the cool- 
ing radius and given k can be used to measure M{r). 

• ROSAT HRI surface brightness data from three Abell 
clusters are best- fit by models in which fc ~ 1. The measured 
core radii are then 110-140 h'^Q kpc. Larger values of k give 
smaller (or non-existent) core radii, but have gas tempera- 
tures which rise sharply at small radii, in contrast with the 
observations. 

It is unfortunate that the core radii deduced for the 
matter distributions for the Abell clusters under investiga- 
tion are only slightly smaller than the radii of the cooling 
fiows. It is possible that the self-similar, steady-state as- 
sumption is a poor one at the edge of the flow (it is also 
possible that it is valid well beyond the cooling radius) and 
that models could be found which were consistent with a 
much larger core radius. 

Ideally, one would like a much better resolution at 
smaller radii, so that the tendency of the emissivity pro- 
file to a constant slope, indicative of the particular value of 
k, could be checked. This should be possible for some of the 
closer clusters, such as Virgo. 
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Note that the sohitions for E and r given in this paper 
do not depend upon the normalisations of the gas and total 
gravitational masses. If desired, these can be determined by 
fixing the overall temperature and luminosity. The analysis 
of Gunn & Thomas (1996) shows that the gas density will 
be slightly lower and the total mass density slightly higher 
than in the equivalent single-phase analysis. 

Cooling flows in individual galaxies are much better 
resolved than those in clusters and so look like promising 
candidates for the kind of modelling discussed here. How- 
ever, the lower temperatures leads to complications such as 
a varying slope, a, for the cooling function, and larger cor- 
rections for absorption and emission lying outside the pass- 
band of the detector. There may also be sources of mass and 
energy injection into the flow. 
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APPENDIX A: DERIVATION OF THE 
MULTIPHASE COOLING FLOW EQUATIONS 

We assume an emulsion of density phases which comove with 
the flow. The distribution is described by the volume frac- 
tion, f{p,r,t), such that f dp is the fractional volume occu- 



pied by phases with densities in the range p to p + dp. Then 
J f dp — 1, and the mean density is p — J fpdp. 
Mass conservation gives 



^(p/) + V.(up/) + ^(pp/) 



0, 



(Al) 



where u is the rate of change of position and p is the rate of 
chan ge o f density following the flow. The final term in equa- 



tion 



Al 



is the equivalent in density space of the divergence 
in velocity space. 

Integrating over all densities we obtain 

^ + \7.u + l3 = 0, (A2) 
P 

where 

P=-lim (ppf). (A3) 

p pi — >oc 

This is equivalent to the usual single-phase equation (e.g. 
Thomas 1988a) except that the mass deposition rate is spec- 
ified in terms of / rather than being a free parameter. 

To find how the volume fraction changes with time we 
use the energy equation, 

ps = -n^A, (A4) 

where s is the entropy and n^A is the radiated power per 
unit volume. For a fully-ionised plasma 



7-1 finiH \ p~' 



(A5) 



where P is the pressure, is the Boltzmann constant, pniH 
is the mass per particle and 7 — 5/3. Then 

p IP 7-1 n^A , . , 

^-7P + V— ^^'^ 
Over a wide temperature range appropriate to clusters 
the cooling function can be approximated by a power-law, 
A oc T", where a ~ 0.5. Then equation A6 can be simplified 



by moving to a new density variable. Writing 
p = po(r,t)w 
we obtain 



(A7) 
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(2- 



po 
po 



IP 



(2- 



,7-lnoA(ro) 



7 



(A8) 



If the adiabatic compression term is removed by setting 
P (X pI, then the energy equation takes a particularly simple 
form, w ^constant. However, a more useful choice is to take 
Po oc p. From the final term of equation Al, we see that at 
high density when cooling is dominant, then ppf ~ constant. 
This motivates the substitution 



(2- 



(4-q)/(2-q) 



Po 



Sr(™,r,t). 



(A9) 



Then, using equations A.1, A6 and we obtain the follow- 
ing equation for the covariant derivate of g (i.e. following 
the fluid flow): 



Q , N Po 2 — a P „ 

^ + 3 - + V.u = 0. 

g po -y P 



(AlO) 
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